In [1]:
# Import Libraries

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import plotly
import plotly.offline as pyoff
import plotly.graph_objs as go
import plotly.express as px
import chart_studio
import chart_studio.plotly as py
import calmap
import datetime
import tensorflow as tf
import os
import random
import re
import plotly.offline as pyoff
import plotly.graph_objs as go
import swifter

from scipy import signal
from datetime import date
from plotly.subplots import make_subplots
from itertools import cycle, product
from statsmodels.tsa.seasonal import STL, seasonal_decompose
from scipy.stats import boxcox
from dateutil.parser import parse
from pmdarima.arima import auto_arima
from pmdarima.utils import diff_inv
from statsmodels.tsa.stattools import adfuller
from sklearn.model_selection import TimeSeriesSplit
from tensorflow.keras.layers import LSTM, Dense, BatchNormalization
from tensorflow.keras import Sequential
from tensorflow.keras.backend import clear_session
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator
from tensorflow.keras.initializers import *
from tensorflow.keras import optimizers
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, mean_absolute_error, r2_score
from sklearn.linear_model import LinearRegression
from scipy.special import boxcox1p, inv_boxcox1p
import matplotlib.patches as mpatches
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.model_selection import GridSearchCV
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from tensorflow.keras.optimizers import Adam
from IPython.display import HTML, display
from swifter import set_defaults
from xgboost import XGBRegressor
from xgboost import plot_importance, plot_tree
In [2]:
# Versões dos pacotes usados neste jupyter notebook
%reload_ext watermark
%watermark -a "Herikc Brecher" --iversions

Author: Herikc Brecher

calmap      : 0.0.9
swifter     : 1.3.4
scipy       : 1.5.4
pandas      : 1.4.2
seaborn     : 0.11.2
matplotlib  : 3.5.1
tensorflow  : 2.10.0
keras       : 2.10.0
chart_studio: 1.1.0
numpy       : 1.20.0
re          : 2.2.1
plotly      : 5.6.0
In [3]:
# Variaveis globais
SEED = 84796315
FEATURES = 7
EPOCHS = 1000
BATCH_SIZE = 1000
EXECUTE_GRID_SEARCH = True
In [4]:
# Configurando seeds
os.environ['PYTHONHASHSEED'] = str(SEED)
tf.random.set_seed(SEED)
np.random.seed(SEED)
random.seed(SEED)
In [5]:
# Exibindo toda tela
display(HTML('<style>.container { width:100% !important; }</style>'))
pd.options.plotting.backend = 'matplotlib'
In [6]:
# Configurando swifter 
set_defaults(
    npartitions = None,
    dask_threshold = 1,
    scheduler = "processes",
    progress_bar = True,
    progress_bar_desc = None,
    allow_dask_on_strings = True,
    force_parallel = True,
)
In [7]:
# Configurando arquivo de API para dar upload aos graficos do plotly
read_file = True

if read_file:
    try:
        with open('api.txt') as f:
            lines = f.readlines()

        # Ler usuario
        username = lines[0].split(':')[1].split('\n')[0]

        # Ler api key
        api_key = lines[1].split(':')[1]
    except:
        read_file = False
In [8]:
# True: Upload dos graficos; False: Não irá dar upload
upload_to_dash_desejado = True

if read_file:
    # Configurando key para upar no dash do plotly
    chart_studio.tools.set_credentials_file(username=username, api_key=api_key)

    # Configurando Privacidade
    chart_studio.tools.set_config_file(world_readable = True, sharing = 'public')

    # Variavel para upar no dash
    upload_to_dash = upload_to_dash_desejado
else:
    upload_to_dash = False

1. Preparando Conjunto de Pedidos¶

In [9]:
# Import dataset
dtOrders = pd.read_csv('../data/olist_orders_dataset.csv', encoding = 'utf8')
In [10]:
# Colunas do tipo data
dateColumns = ['order_purchase_timestamp', 'order_approved_at', 'order_delivered_carrier_date',\
               'order_delivered_customer_date', 'order_estimated_delivery_date']

# Dataset de analise temporal
dtOrdersAdjusted = dtOrders.copy()
In [11]:
# Convertendo columas de data para date
for col in dateColumns:
    dtOrdersAdjusted[col] = pd.to_datetime(dtOrdersAdjusted[col], format = '%Y-%m-%d %H:%M:%S')
In [12]:
# Dropando valores NA
dtOrdersAdjusted = dtOrdersAdjusted.dropna()
In [13]:
dtOrdersAdjusted.dtypes
Out[13]:
order_id                                 object
customer_id                              object
order_status                             object
order_purchase_timestamp         datetime64[ns]
order_approved_at                datetime64[ns]
order_delivered_carrier_date     datetime64[ns]
order_delivered_customer_date    datetime64[ns]
order_estimated_delivery_date    datetime64[ns]
dtype: object

2. Iniciando Analise Seasonal¶

In [14]:
dtHistory = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date

start = dtHistory.min()
end = dtHistory.max()

idx = pd.date_range(start, end, normalize = True)

seriesOriginal = dtHistory.value_counts(sort = False).sort_index().reindex(idx, fill_value = 0)

dtHistory = pd.DataFrame(seriesOriginal).reset_index()

Principais outliers identificados:

  • 1 de setembro de 2016 a 31 de dezembro de 2016: Dados quase inexistentes
  • 24 de novembro de 2017: Pico de venda devido ao evento da blackfriday
  • 17 de agosto de 2017 a 17 de outubro de 2017: Queda repentina nos dados
In [15]:
dtHistory.rename(columns = {'index': 'Data', 'order_purchase_timestamp': 'Vendas'}, inplace = True)
In [16]:
dtHistory
Out[16]:
Data Vendas
0 2016-09-15 1
1 2016-09-16 0
2 2016-09-17 0
3 2016-09-18 0
4 2016-09-19 0
... ... ...
709 2018-08-25 69
710 2018-08-26 73
711 2018-08-27 66
712 2018-08-28 39
713 2018-08-29 11

714 rows × 2 columns

In [17]:
# Plot

# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['Data'],
                        y = dtHistory['Vendas'])]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
                        yaxis = {'title': 'Vendas'},
                        title = 'Vendas por dia')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.update_layout(yaxis_range = [0, 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'vendas_dia', auto_open = False)
In [18]:
# Remove outliers
seriesOriginal = seriesOriginal[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]
pred_range = pd.date_range(datetime.date(2018, 8, 17), datetime.date(2018, 10, 17))
dtHistory = pd.DataFrame(seriesOriginal).reset_index()
In [19]:
dtHistory.rename(columns = {'index': 'Data', 'order_purchase_timestamp': 'Vendas'}, inplace = True)
In [20]:
# Plot

# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['Data'],
                        y = dtHistory['Vendas'])]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
                        yaxis = {'title': 'Vendas'},
                        title = 'Vendas por Dia')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.update_layout(yaxis_range = [0, 800])

pyoff.iplot(fig) 

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'vendas_dia_sem_outlier', auto_open = False)
In [21]:
# Plot

# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['Data'],
                        y = dtHistory['Vendas'],
                        showlegend = False),
            go.Scatter(x = dtHistory['Data'],
                       y = -dtHistory['Vendas'],
                       fill='tonexty',
                       showlegend = False)]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
                        yaxis = {'title': 'Vendas'},
                        title = 'Vendas por Dia Expandida')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.update_traces(line_color = 'green')
fig.update_layout(yaxis_range = [-800 , 800])

pyoff.iplot(fig) 

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'vendas_dia_expandida', auto_open = False)
In [22]:
#Plot histórico de vendas por dia, mês e ano
fig, caxs = calmap.calendarplot(seriesOriginal, daylabels = 'MTWTFSS', fillcolor = 'grey',cmap = 'YlGn', fig_kws = dict(figsize = (18, 9)))
fig.suptitle('Histórico de Vendas', fontsize = 22)

fig.subplots_adjust(right = 0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.03, 0.67])
fig.colorbar(caxs[0].get_children()[1], cax = cbar_ax)

plt.show()
In [23]:
# Criar grafico na estrutura STL 4 layers
def add_stl_plot(fig, res, legend):
    axs = fig.get_axes()
    
    # Nome de cada um dos subplots
    comps = ['trend', 'seasonal', 'resid']
    for ax, comp in zip(axs[1:], comps):
        series = getattr(res, comp)
        if comp == 'resid':
            ax.plot(series, marker = 'o', linestyle = 'none')
        else:
            ax.plot(series)
            ax.legend(legend, frameon = False)
In [24]:
# Gerar STL
stl = STL(seriesOriginal)
stl_res = stl.fit()
fig = stl_res.plot()
fig.set_size_inches((20, 12))
plt.show()
In [25]:
# Gerar STL não robusto e concatenar ao robusto
stl = STL(seriesOriginal, robust = True)
res_robust = stl.fit()
fig = res_robust.plot()
fig.set_size_inches((20, 12))
res_non_robust = STL(seriesOriginal, robust = False).fit()
add_stl_plot(fig, res_non_robust, ['Robusto', 'Não Robusto'])
In [26]:
# Additive Decomposition
additive_decomposition = seasonal_decompose(seriesOriginal, model = 'additive', period = FEATURES)

# Plot
additive_decomposition.plot()
plt.tight_layout(rect = [0, 0.03, 1, 0.95])
plt.rcParams["figure.figsize"] = (20, 10)

plt.show()
In [27]:
# Removendo tendência
detrended = signal.detrend(seriesOriginal)

# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
                        y = detrended)]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'Vendas Diárias Removendo a Tendência')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.update_layout(yaxis_range = [0 , 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'vendas_diarias_sem_trend', auto_open = False)
In [28]:
# Gerando STL para separar cada um dos componentes
stl = STL(seriesOriginal)
res = stl.fit()

# Separando seriesDeseasonal
seriesDeseasonal = res.observed - res.seasonal

# Separando boxcox
seriesBoxCox, lmbda = boxcox(seriesOriginal + 1)
seriesBoxCox = pd.Series(seriesBoxCox, index = seriesOriginal.index)

# Separando stationary
seriesResidual = seriesOriginal.diff(7).dropna()

2.1 Teste Estacionário ADF¶

Os testes abaixo concluiram:

O teste aceita a hipótese nula em que a série não é estácionária para os dados originais e deseasonal. Já para os dados residuais esses aceitaram a hipótese alternativa que os dados são estacionários.

ADF teste:

  • Hipótese Nula(HO): A série possui unit root ou não é estacionária.
  • Hipótese Alternativa(HA): A série não possui unit root ou é estacionária.
In [29]:
print("Os dados são estacionários?\n")
testResult = adfuller(seriesOriginal, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")

for p, v in testResult[4].items():
    print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))

Os dados são estacionários?

Valor Teste = -2.616
Valor de P: = 0.090

Valores Críticos:
	1%: -3.441694608475642 - O dataset não é estacionário com 99% de confiança
	5%: -2.866544718556839 - O dataset não é estacionário com 95% de confiança
	10%: -2.5694353738653684 - O dataset  é estacionário com 90% de confiança
In [30]:
print("Os dados deseasonal são estacionários?")
testResult = adfuller(seriesDeseasonal, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")

for p, v in testResult[4].items():
    print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))

Os dados deseasonal são estacionários?
Valor Teste = -2.536
Valor de P: = 0.107

Valores Críticos:
	1%: -3.441694608475642 - O dataset não é estacionário com 99% de confiança
	5%: -2.866544718556839 - O dataset não é estacionário com 95% de confiança
	10%: -2.5694353738653684 - O dataset não é estacionário com 90% de confiança
In [31]:
print("Os dados boxcox são estacionários?")
testResult = adfuller(seriesBoxCox, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")

for p, v in testResult[4].items():
    print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))

Os dados boxcox são estacionários?
Valor Teste = -2.874
Valor de P: = 0.048

Valores Críticos:
	1%: -3.441694608475642 - O dataset não é estacionário com 99% de confiança
	5%: -2.866544718556839 - O dataset  é estacionário com 95% de confiança
	10%: -2.5694353738653684 - O dataset  é estacionário com 90% de confiança
In [32]:
print("Os dados residuais são estacionários?")
testResult = adfuller(seriesResidual, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")

for p, v in testResult[4].items():
    print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))

Os dados residuais são estacionários?
Valor Teste = -6.802
Valor de P: = 0.000

Valores Críticos:
	1%: -3.441834071558759 - O dataset  é estacionário com 99% de confiança
	5%: -2.8666061267054626 - O dataset  é estacionário com 95% de confiança
	10%: -2.569468095872659 - O dataset  é estacionário com 90% de confiança

3. Modelagem¶

Toda a etapa de modelagem será considerada com 5 passos a frente de previsão.

In [33]:
# Controle de resultados de toda fase de modelagem
result = pd.DataFrame(columns = ['Algorithm', 'MSE', 'RMSE', 'MAE', 'MAPE', 'Mean_True_Value', 'Mean_Predict_Value'])
In [34]:
split_range = TimeSeriesSplit(n_splits = 8, max_train_size = pred_range.shape[0], test_size = pred_range.shape[0])
In [35]:
# Adiciona o registro ao dataset
def record(result, algorithm, mse = -1, rmse = -1, mae = -1, mape = -1, mrv = -1, mpv = -1, show = True):
    new = pd.DataFrame(dict(Algorithm = algorithm, MSE = mse, RMSE = rmse, MAE = mae, MAPE = mape, Mean_True_Value = mrv,\
                            Mean_Predict_Value = mpv), index = [0])
    result = pd.concat([result, new], ignore_index = True)
    
    if show:
        display(result)
    
    return result
In [36]:
# Plot no formato de 4 layers, seguindo o STL para cada um dos modelos
def plot(index, pred, mse, title, fig = None, ax = None, ylabel = ''):
    global seriesOriginal
    
    empty_fig = fig is None
    
    if empty_fig:
        fig, ax = plt.subplots(figsize = (13, 6))
    else: 
        ax.set_ylabel(ylabel)
                
    ax.set_title(title)    
    patch_ = mpatches.Patch(color = 'white', label = f'MSE: {np.mean(mse):.1e}')
    L1 = ax.legend(handles = [patch_], loc = 'upper left', fancybox = True, framealpha = 0.7,  handlelength = 0)
    ax.add_artist(L1)
    
    sns.lineplot(x = seriesOriginal.index, y = seriesOriginal, label = 'Real', ax = ax)
    sns.lineplot(x = index, y = pred, label = 'Previsto', ax = ax)
    ax.axvline(x = index[0], color = 'red')
    ax.legend(loc = 'upper right')
    ax.set_ylim([0, 800])
    
    if empty_fig:
        plt.show()
    else:
        return fig
In [37]:
# Calculo para previsão e teste quando utilizado a série Original
def calcPredTestOriginal(train, pred, test):
    return pred, test, 0
In [38]:
# Calculo para previsão e teste quando utilizado a série seriesDeseasonal
def calcPredTestseriesDeseasonal(train, pred, test):
    # Removendo a sazonalidade da série e convertendo para o shape correto
    last_seasonal = res.seasonal.reindex_like(train).tail(stl.period)
    pred = pred + np.fromiter(cycle(last_seasonal), count = pred.shape[0], dtype = float)
    test = test + res.seasonal.reindex_like(test)
    
    return pred, test, 1
In [39]:
# Calculo para previsão e teste quando utilizado a série BoxCox
def calcPredTestBoxCox(train, pred, test):
    # Reverdendo a normalização do boxcox
    pred = inv_boxcox1p(pred, lmbda)
    test = inv_boxcox1p(test, lmbda)
    
    return pred, test, 2
In [40]:
# Calculo para previsão e teste quando utilizado a série Stationary
def calcPredTestStationary(train, pred, test):
    # Calculando a diferença da sazonalidade
    xi = seriesOriginal.reindex_like(train).tail(FEATURES)
    
    totalLen = len(pred) + len(xi) 
    ix = pd.date_range(xi.index[0], periods = totalLen)  
    inv = diff_inv(pred, FEATURES, xi = xi) + np.fromiter(cycle(xi), count = totalLen, dtype = float)  
    inv = pd.Series(inv, index = ix, name = 'Vendas')
    pred = inv.iloc[FEATURES:]
    
    totalLen = len(test) + len(xi) 
    ix = pd.date_range(xi.index[0], periods = totalLen)  
    inv = diff_inv(test, FEATURES, xi = xi) + np.fromiter(cycle(xi), count = totalLen, dtype = float)  
    inv = pd.Series(inv, index = ix, name = 'Vendas')
    test = inv.iloc[FEATURES:]
    
    return pred, test, 3
In [41]:
def calculate_metrics_off(test, pred):
    mse = mean_squared_error(test, pred, squared = True)
    rmse = mean_squared_error(test, pred, squared = False)
    mae = mean_absolute_error(test, pred)
    mape = mean_absolute_percentage_error(test, pred) * 100
    mtv = np.mean(test)
    mpv = np.mean(pred)
    
    return mse, rmse, mae, mape, mtv, mpv
In [42]:
def calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv):
    mse_, rmse_, mae_, mape_, mtv_, mpv_ = calculate_metrics_off(test, pred)
    
    mse.append(mse_)
    rmse.append(rmse_)
    mae.append(mae_)
    mape.append(mape_)
    mtv.append(mtv_)
    mpv.append(mpv_)
    
    return mse, rmse, mae, mape, mtv, mpv

3.1 TSR (Time Series Regression)¶

In [43]:
# Report para Time Series Regressor, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportTSR(data, modelName, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mape = []
    mtv = []
    mpv = []
    
    title = modelName + ' - Time Series Regression'
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
    
        gen = TimeseriesGenerator(train, train, FEATURES, batch_size = BATCH_SIZE)

        X_train = gen[0][0]
        y_train = gen[0][1]

        lr = LinearRegression()
        lr.fit(X_train, y_train)
        X_pred = y_train[-FEATURES:].reshape(1,-1)
        pred = np.empty(test.shape[0])

        for i in range(len(pred)):
            forecast = lr.predict(X_pred)
            X_pred = np.delete(X_pred, 0, 1)
            X_pred = np.concatenate((X_pred, forecast.reshape(-1, 1)), 1)    
            pred[i] = forecast
        
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
In [44]:
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.tight_layout()
plt.close()
In [45]:
reportTSR(seriesOriginal.copy(), 'Original', calcPredTestOriginal)
Out[45]:

3.2 Deseasonal - TSR (Time Series Regression)¶

In [46]:
reportTSR(seriesDeseasonal.copy(), 'Deseasonal', calcPredTestseriesDeseasonal)
Out[46]:

3.3 BoxCox - TSR (Time Series Regression)¶

In [47]:
reportTSR(seriesBoxCox.copy(), 'BoxCox', calcPredTestBoxCox)
Out[47]:

3.4 Residual - TSR (Time Series Regression)¶

In [48]:
reportTSR(seriesResidual.copy(), 'Stationary', calcPredTestStationary)
Out[48]:
In [49]:
result
Out[49]:
Algorithm MSE RMSE MAE MAPE Mean_True_Value Mean_Predict_Value
0 Original - Time Series Regression 5296.329792 60.784585 42.960631 25.553663 179.560484 168.008596
1 Deseasonal - Time Series Regression 5081.77776 58.208746 40.122247 23.488563 179.560484 164.670085
2 BoxCox - Time Series Regression 5328.602465 61.060147 43.254815 25.416806 179.560484 165.399796
3 Stationary - Time Series Regression 6057.828087 64.323443 46.538896 27.899678 179.560484 184.267746

3.5 Exponential Smoothing¶

In [50]:
# Função utilizada para o hypertuning de alpha, beta e gamma do Exponential Smoothing
def GSES(data, modelName, alpha, beta, gamma, damping_trend, calcFunction):    
    mse = []
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        try:
            with catch_warnings():
                filterwarnings('ignore')
                ES = (
                    ExponentialSmoothing(train, trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
                    .fit(smoothing_level = alpha, smoothing_trend = beta, smoothing_seasonal = gamma, method = 'ls', damping_trend = damping_trend)
                )

                pred = ES.forecast(test.shape[0])

                pred, test, _ = calcFunction(train, pred, test)

                mse.append(mean_squared_error(pred, test, squared = True))
        
        except:
            mse.append(-1)
    
    return np.mean(mse)
In [51]:
# Função utilizada para o hypertuning de demais parâmetros do Exponential Smoothing
def GSESOPT(data, modelName, trend, season, periods, bias, method, calcFunction):
    mse = []
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        try:
            with catch_warnings():
                filterwarnings('ignore')
                ES = (
                    ExponentialSmoothing(train, trend = trend, seasonal = season, seasonal_periods = periods)
                    .fit(remove_bias = bias, method = method, optimized = True)
                )

                pred = ES.forecast(test.shape[0])

                pred, test, _ = calcFunction(train, pred, test)

                mse.append(mean_squared_error(pred, test, squared = True))       
        except:
            mse.append(-1)
    
    return np.mean(mse)
In [52]:
# Report para Exponential Smoothing, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportES(data, modelName, model_kwargs, fit_kwargs, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mape = []
    mtv = []
    mpv = []
    
    title = modelName + ' - Exponential Smoothing'
    indexPlot = 0
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        ES = (
            ExponentialSmoothing(train, trend = model_kwargs['trend'], seasonal = model_kwargs['seasonal'], seasonal_periods = FEATURES, damped_trend = model_kwargs['damped_trend'])
            .fit(smoothing_level = fit_kwargs['smoothing_level'], smoothing_trend = fit_kwargs['smoothing_trend'],\
                 smoothing_seasonal = fit_kwargs['smoothing_seasonal'], method = fit_kwargs['method'], damping_trend = fit_kwargs['damping_trend'])
        )
        
        pred = ES.forecast(test.shape[0])
             
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
In [53]:
# Função para gerar tabela de hypertuning ampla
def exp_smoothing_configs(seasonal = [None]):
    models = list()
    # Lista de argumentos
    t_params = ['add', 'mul']
    s_params = ['add', 'mul']
    p_params = seasonal
    r_params = [True, False]
    method_params = ['L-BFGS-B' , 'TNC', 'SLSQP', 'Powell', 'trust-constr', 'bh', 'ls']
    
    # Gerando lista de argumentos
    for t in t_params:
        for s in s_params:
            for p in p_params:
                for r in r_params:
                    for m in method_params:
                        cfg = [t, s, p, r, m]
                        models.append(cfg)
    return models
In [54]:
%%time

# Treinamento do modelo 
if EXECUTE_GRID_SEARCH:
    # Gerando tabela de hypertunning
    alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
    hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesOriginal.copy(), 'Original',\
                                                x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestOriginal), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Dask Apply:   0%|          | 0/48 [00:00<?, ?it/s]
alpha beta gamma damping_trend mse
1808 0.1 0.8 0.0 0.8 4749.165741 
CPU times: total: 26.7 s
Wall time: 1min 55s
In [55]:
%%time

# Se True irá treinar com a nova lista mais ampla (pode demorar)
if EXECUTE_GRID_SEARCH:
    # Criando lista de argumentos ampla
    params_ = exp_smoothing_configs([FEATURES])
    
    # Convertendo dicionário de argumentos para dataframe
    hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
    
    # Executando hypertuning
    hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesOriginal.copy(), 'Original',\
                                                     x.trend, x.season, x.periods, x.bias, x.method, calcPredTestOriginal),\
                                   axis = 1)
    
    # Busca em query pelos melhores argumentos
    display(hyperparam_.query('mse == mse.min() and mse != -1'))

Dask Apply:   0%|          | 0/28 [00:00<?, ?it/s]
trend season periods bias method mse
16 add mul 7 True SLSQP 5217.975062 
CPU times: total: 2.42 s
Wall time: 1min 8s
In [56]:
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13, 6))
figs.align_ylabels()
figs.tight_layout()
plt.close()
In [57]:
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.1, smoothing_trend = 0.8, smoothing_seasonal = 0, method = 'ls', damping_trend = 0.8)
In [58]:
reportES(seriesOriginal.copy(), 'Original', model_kwargs, fit_kwargs, calcPredTestOriginal)
Out[58]:

3.6 Deseasonal - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesOriginal para seriesDeseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [59]:
%%time

if EXECUTE_GRID_SEARCH:
    alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
    hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesDeseasonal.copy(), 'seriesDeseasonal',\
                                                x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestseriesDeseasonal), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Dask Apply:   0%|          | 0/48 [00:00<?, ?it/s]
alpha beta gamma damping_trend mse
1258 0.1 0.2 0.5 0.8 4774.382729 
CPU times: total: 27.7 s
Wall time: 1min 57s
In [60]:
%%time
if EXECUTE_GRID_SEARCH:
    params_ = exp_smoothing_configs([FEATURES])
    hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
    
    hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesDeseasonal.copy(), 'seriesDeseasonal',\
                                                     x.trend, x.season, x.periods, x.bias, x.method, calcPredTestseriesDeseasonal),\
                                   axis = 1)
    
    display(hyperparam_.query('mse == mse.min() and mse != -1'))

Dask Apply:   0%|          | 0/28 [00:00<?, ?it/s]
trend season periods bias method mse
34 mul add 7 True ls 5179.581951 
CPU times: total: 2.64 s
Wall time: 1min 51s
In [61]:
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.1, smoothing_trend = 0.2, smoothing_seasonal = 0.5, method = 'ls', damping_trend = 0.8)
In [62]:
reportES(seriesDeseasonal.copy(), 'Deseasonal', model_kwargs, fit_kwargs, calcPredTestseriesDeseasonal)
Out[62]:

3.7 BoxCox - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesOriginal para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [63]:
%%time

if EXECUTE_GRID_SEARCH:
    alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
    hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesBoxCox.copy(), 'BoxCox',\
                                                x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestBoxCox), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Pandas Apply:   0%|          | 0/10000 [00:00<?, ?it/s]
alpha beta gamma damping_trend mse
1708 0.1 0.7 0.0 0.8 4653.246559 
CPU times: total: 11min 33s
Wall time: 11min 37s
In [64]:
%%time
if EXECUTE_GRID_SEARCH:
    params_ = exp_smoothing_configs([FEATURES])
    hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
    
    hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesBoxCox.copy(), 'BoxCox',\
                                                  x.trend, x.season, x.periods, x.bias, x.method, calcPredTestBoxCox),\
                                   axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Pandas Apply:   0%|          | 0/56 [00:00<?, ?it/s]
alpha beta gamma damping_trend mse
1708 0.1 0.7 0.0 0.8 4653.246559 
CPU times: total: 3min 59s
Wall time: 4min 4s
In [65]:
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.1, smoothing_trend = 0.7, smoothing_seasonal = 0.0, method = 'ls', damping_trend = 0.8)
In [66]:
reportES(seriesBoxCox.copy(), 'BoxCox', model_kwargs, fit_kwargs, calcPredTestBoxCox)
Out[66]:

3.8 Residual - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesOriginal para seriesResidual, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [67]:
%%time
if EXECUTE_GRID_SEARCH:
    alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
    hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesResidual.copy(), 'Stationary',\
                                                x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestStationary), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Dask Apply:   0%|          | 0/48 [00:00<?, ?it/s]
alpha beta gamma damping_trend mse
712 0.0 0.7 0.1 0.2 5658.627132 
CPU times: total: 31.6 s
Wall time: 2min 13s
In [68]:
%%time
if EXECUTE_GRID_SEARCH:
    params_ = exp_smoothing_configs([FEATURES])
    
    hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
    
    hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesResidual.copy(), 'Stationary',\
                                                  x.trend, x.season, x.periods, x.bias, x.method, calcPredTestStationary),\
                                   axis = 1)
    
    display(hyperparam_.query('mse == mse.min() and mse != -1'))

Dask Apply:   0%|          | 0/28 [00:00<?, ?it/s]
trend season periods bias method mse 
CPU times: total: 172 ms
Wall time: 1min 24s
In [69]:
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.0, smoothing_trend = 0.2, smoothing_seasonal = 0.1, method = 'ls', damping_trend = 0.2)
In [70]:
reportES(seriesResidual.copy(), 'Stationary', model_kwargs, fit_kwargs, calcPredTestStationary)
Out[70]:
In [71]:
result
Out[71]:
Algorithm MSE RMSE MAE MAPE Mean_True_Value Mean_Predict_Value
0 Original - Time Series Regression 5296.329792 60.784585 42.960631 25.553663 179.560484 168.008596
1 Deseasonal - Time Series Regression 5081.77776 58.208746 40.122247 23.488563 179.560484 164.670085
2 BoxCox - Time Series Regression 5328.602465 61.060147 43.254815 25.416806 179.560484 165.399796
3 Stationary - Time Series Regression 6057.828087 64.323443 46.538896 27.899678 179.560484 184.267746
4 Original - Exponential Smoothing 4749.165741 55.864659 38.516078 22.095183 179.560484 165.919134
5 Deseasonal - Exponential Smoothing 4774.382729 55.5275 37.488666 21.8423 179.560484 166.749778
6 BoxCox - Exponential Smoothing 4653.246555 54.946508 37.555881 21.63599 179.560484 166.093056
7 Stationary - Exponential Smoothing 5658.627363 63.86568 43.658278 26.080376 179.560484 179.444888

3.9 ARIMA¶

In [72]:
# Report do algoritmo arima, também é adicionado a base de resultados e realizado o plot de acompanhamento
def reportArima(arimaModel, modelName, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mape = []
    mtv = []
    mpv = []
    
    title = modelName + ' - Auto Arima'
    indexPlot = 0
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        arimaModel.fit(train)
        pred = arimaModel.predict(test.shape[0])
             
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
In [73]:
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()
In [74]:
# Correlação entre os periodos com ARIMA

lags = 90

with catch_warnings():
    filterwarnings('ignore')
    fig, ax = plt.subplots(2, figsize = (12, 6), sharex = True)
    plot_acf(seriesOriginal.diff().dropna(), ax = ax[0], lags = lags, missing = 'drop')
    plot_pacf(seriesOriginal.diff().dropna(), ax = ax[1], lags = lags)
    plt.show()
In [75]:
%%time

# Utilizando o auto arima para descobrir os argumentos ideias baseados no conjunto de dado informado
data = seriesOriginal.copy()
arimaModel = auto_arima(seriesOriginal.copy(), m = FEATURES, seasonal = False)
arimaModel

C:\Users\herik\anaconda3\lib\site-packages\pmdarima\arima\_validation.py:62: UserWarning:

m (7) set for non-seasonal fit. Setting to 0


CPU times: total: 2.28 s
Wall time: 2.3 s
Out[75]:
ARIMA(order=(1, 1, 1), scoring_args={}, suppress_warnings=True,
      with_intercept=False)
In [76]:
reportArima(arimaModel, 'Original', calcPredTestOriginal)
Out[76]:

3.10 Deseasonal - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada de seriesOriginal para seriesDeseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [77]:
%%time
data = seriesDeseasonal.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = False)
arimaModel

C:\Users\herik\anaconda3\lib\site-packages\pmdarima\arima\_validation.py:62: UserWarning:

m (7) set for non-seasonal fit. Setting to 0


CPU times: total: 2.59 s
Wall time: 2.61 s
Out[77]:
ARIMA(order=(2, 1, 1), scoring_args={}, suppress_warnings=True)
In [78]:
reportArima(arimaModel, 'Deseasonal', calcPredTestseriesDeseasonal)
Out[78]:

3.11 BoxCox - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada de seriesOriginal para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [79]:
%%time
data = seriesBoxCox.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = True)
arimaModel

CPU times: total: 26.5 s
Wall time: 26.6 s
Out[79]:
ARIMA(order=(1, 1, 1), scoring_args={}, seasonal_order=(1, 0, 1, 7),
      suppress_warnings=True)
In [80]:
reportArima(arimaModel, 'BoxCox', calcPredTestBoxCox)
Out[80]:

3.12 Residual - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada original para stationary, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [81]:
%%time
data = seriesResidual.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = True)
arimaModel

CPU times: total: 43.6 s
Wall time: 43.7 s
Out[81]:
ARIMA(order=(1, 0, 5), scoring_args={}, seasonal_order=(0, 0, 1, 7),
      suppress_warnings=True)
In [82]:
reportArima(arimaModel, 'Stationary', calcPredTestStationary)
Out[82]:
In [83]:
result
Out[83]:
Algorithm MSE RMSE MAE MAPE Mean_True_Value Mean_Predict_Value
0 Original - Time Series Regression 5296.329792 60.784585 42.960631 25.553663 179.560484 168.008596
1 Deseasonal - Time Series Regression 5081.77776 58.208746 40.122247 23.488563 179.560484 164.670085
2 BoxCox - Time Series Regression 5328.602465 61.060147 43.254815 25.416806 179.560484 165.399796
3 Stationary - Time Series Regression 6057.828087 64.323443 46.538896 27.899678 179.560484 184.267746
4 Original - Exponential Smoothing 4749.165741 55.864659 38.516078 22.095183 179.560484 165.919134
5 Deseasonal - Exponential Smoothing 4774.382729 55.5275 37.488666 21.8423 179.560484 166.749778
6 BoxCox - Exponential Smoothing 4653.246555 54.946508 37.555881 21.63599 179.560484 166.093056
7 Stationary - Exponential Smoothing 5658.627363 63.86568 43.658278 26.080376 179.560484 179.444888
8 Original - Auto Arima 5707.329361 63.881702 45.981605 27.441516 179.560484 169.582966
9 Deseasonal - Auto Arima 9591.527185 77.052424 59.641758 35.25454 179.560484 191.567949
10 BoxCox - Auto Arima 11358.34545 84.301119 65.970324 39.351481 179.560484 206.5619
11 Stationary - Auto Arima 9014.243165 73.755288 56.443311 33.402952 179.560484 192.709473

3.13 LSTM¶

In [84]:
# Redefinindo variaveis globais para o treino utilizando LSTM

BATCH_SIZE = 30
In [85]:
# hypertuning do algoritmo de LSTM
def GSLSTM(data, learning_rate, calcFunction):
    mse = []
    
    # Crossvalidation para cada parte do conjunto
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]

        try:
            with catch_warnings():
                filterwarnings('ignore')
                
                # Normalização e reshape do conjunto de treino
                ss = StandardScaler()
                ss.fit(train.values.reshape(-1, 1))
                train_input = ss.transform(train.values.reshape(-1, 1))
                
                # Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
                test_input = train_input[-(FEATURES + 1):]
                test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
                train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)
                
                # Reset da sessão
                clear_session()
                
                # Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
                initializer = GlorotUniform(seed = SEED)
                model = Sequential()
                
                # 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
                model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
                model.add(Dense(1, kernel_initializer = initializer))
                model.add(BatchNormalization())
                model.add(Dense(1, kernel_initializer = initializer))
                
                # Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
                early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
                
                # Treinando o modelo com otimizador Adam
                model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
                h = model.fit(train_gen, epochs = EPOCHS, callbacks = [early_stopping], verbose = False)
                pred = np.empty(test.shape[0])
                
                # Realizando predições no conjunto de teste
                for i in range(len(pred)):
                    prediction = model.predict(test_gen, verbose = False)
                    pred[i] = prediction
                    test_input = np.delete(test_input, 0, 0)
                    test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
                    test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
                
                # Reorganizando o shape e chamando a função de calculo
                pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
                pred, test, _ = calcFunction(train, pred, test)      

                mse.append(mean_squared_error(pred, test))
                
        except:
            mse.append(-1)
        
    return np.mean(mse)
In [86]:
# Report do algoritmo LSTM
def reportLSTM(data, modelName, calcFunction, learning_rate):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mape = []
    mtv = []
    mpv = []
    
    title = modelName + ' - Long Short Term Memory (LSTM)'
    
    # Crossvalidation para cada parte do conjunto
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
    
        # Normalização e reshape do conjunto de treino
        ss = StandardScaler()
        ss.fit(train.values.reshape(-1, 1))
        train_input = ss.transform(train.values.reshape(-1, 1))

        # Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
        test_input = train_input[-(FEATURES + 1):]
        test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
        train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)

        # Reset da sessão
        clear_session()
        
        # Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
        initializer = GlorotUniform(seed = SEED)
        model = Sequential()
        
        # 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
        model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
        model.add(Dense(1, kernel_initializer = initializer))
        model.add(BatchNormalization())
        model.add(Dense(1, kernel_initializer = initializer))
        
        # Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
        early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
        
        # Treinando o modelo com otimizador Adam
        model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
        h = model.fit(train_gen, epochs = EPOCHS, callbacks = [early_stopping], verbose = False)
        pred = np.empty(test.shape[0])

        # Realizando predições no conjunto de teste
        for i in range(len(pred)):
            prediction = model.predict(test_gen, verbose = False)
            pred[i] = prediction
            test_input = np.delete(test_input, 0, 0)
            test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
            test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)

        # Reorganizando o shape e chamando a função de calculo
        pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
In [87]:
%%time
if EXECUTE_GRID_SEARCH:
    # Gerando tabela de hypertunning com taxas de learning_rate
    learning_rates = np.logspace(-5, 1, 7)
    hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesOriginal.copy(), x.learning_rate, calcPredTestOriginal), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Pandas Apply:   0%|          | 0/7 [00:00<?, ?it/s]
learning_rate mse
0 0.00001 5506.010033 
CPU times: total: 18min 46s
Wall time: 13min 43s
In [88]:
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()
In [89]:
reportLSTM(seriesOriginal.copy(), 'Original', calcPredTestOriginal, 0.0001)
Out[89]:

3.14 Deseasonal - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para seriesDeseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [90]:
%%time

if EXECUTE_GRID_SEARCH:
    learning_rates = np.logspace(-5, 1, 7)
    hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesDeseasonal.copy(), x.learning_rate, calcPredTestseriesDeseasonal), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Pandas Apply:   0%|          | 0/7 [00:00<?, ?it/s]
learning_rate mse
3 0.01 4868.290229 
CPU times: total: 17min 6s
Wall time: 12min 48s
In [91]:
reportLSTM(seriesDeseasonal.copy(), 'Deseasonal', calcPredTestseriesDeseasonal, 0.01)
Out[91]:

3.15 BoxCox - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [92]:
%%time

if EXECUTE_GRID_SEARCH:
    learning_rates = np.logspace(-5, 1, 7)
    hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesBoxCox.copy(), x.learning_rate, calcPredTestBoxCox), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Pandas Apply:   0%|          | 0/7 [00:00<?, ?it/s]
learning_rate mse
2 0.001 5492.933495 
CPU times: total: 19min 10s
Wall time: 13min 20s
In [93]:
reportLSTM(seriesBoxCox.copy(), 'BoxCox', calcPredTestBoxCox, 0.001)
Out[93]:

3.16 Residual - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para stationary, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

In [94]:
%%time

if EXECUTE_GRID_SEARCH:
    learning_rates = np.logspace(-5, 1, 7)
    hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
    
    hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesResidual.copy(), x.learning_rate, calcPredTestStationary), axis = 1)
    
    display(hyperparam.query('mse == mse.min() and mse != -1'))

Pandas Apply:   0%|          | 0/7 [00:00<?, ?it/s]
learning_rate mse
2 0.001 7389.775821 
CPU times: total: 16min 21s
Wall time: 12min 32s
In [95]:
reportLSTM(seriesResidual.copy(), 'Stationary', calcPredTestStationary, 0.00001)
Out[95]:
In [96]:
result
Out[96]:
Algorithm MSE RMSE MAE MAPE Mean_True_Value Mean_Predict_Value
0 Original - Time Series Regression 5296.329792 60.784585 42.960631 25.553663 179.560484 168.008596
1 Deseasonal - Time Series Regression 5081.77776 58.208746 40.122247 23.488563 179.560484 164.670085
2 BoxCox - Time Series Regression 5328.602465 61.060147 43.254815 25.416806 179.560484 165.399796
3 Stationary - Time Series Regression 6057.828087 64.323443 46.538896 27.899678 179.560484 184.267746
4 Original - Exponential Smoothing 4749.165741 55.864659 38.516078 22.095183 179.560484 165.919134
5 Deseasonal - Exponential Smoothing 4774.382729 55.5275 37.488666 21.8423 179.560484 166.749778
6 BoxCox - Exponential Smoothing 4653.246555 54.946508 37.555881 21.63599 179.560484 166.093056
7 Stationary - Exponential Smoothing 5658.627363 63.86568 43.658278 26.080376 179.560484 179.444888
8 Original - Auto Arima 5707.329361 63.881702 45.981605 27.441516 179.560484 169.582966
9 Deseasonal - Auto Arima 9591.527185 77.052424 59.641758 35.25454 179.560484 191.567949
10 BoxCox - Auto Arima 11358.34545 84.301119 65.970324 39.351481 179.560484 206.5619
11 Stationary - Auto Arima 9014.243165 73.755288 56.443311 33.402952 179.560484 192.709473
12 Original - Long Short Term Memory (LSTM) 48171.637513 139.146814 101.418725 70.382288 179.560484 137.450098
13 Deseasonal - Long Short Term Memory (LSTM) 4989.531811 56.707357 38.964999 22.768957 179.560484 165.563544
14 BoxCox - Long Short Term Memory (LSTM) 5341.400913 62.574598 44.292526 26.569442 179.560484 171.38447
15 Stationary - Long Short Term Memory (LSTM) 16091.513901 94.546789 74.40401 43.338943 179.560484 220.839974

3.16 Original - XGBoost Regressor¶

In [97]:
def create_features(df, label = None):
    # Separando variaveis temporais
    df['date'] = df.index
    df['hour'] = df['date'].dt.hour
    df['dayofweek'] = df['date'].dt.dayofweek
    df['quarter'] = df['date'].dt.quarter
    df['month'] = df['date'].dt.month
    df['year'] = df['date'].dt.year
    df['dayofyear'] = df['date'].dt.dayofyear
    df['dayofmonth'] = df['date'].dt.day
    df['weekofyear'] = df['date'].dt.weekofyear
    
    X = df[['hour','dayofweek','quarter','month','year',
           'dayofyear','dayofmonth','weekofyear']]
    if label:
        y = df[label]
        return X, y
    return X
In [98]:
def reportModeloRegressao(modelo, x_teste, y_teste):  
    y_pred = modelo.predict(x_teste)
    
    erros = abs(y_pred - y_teste)
    mape = 100 * np.mean(erros / y_teste)
    acuracia = 100 - mape
    mse = mean_squared_error(y_teste, y_pred, squared = True)
    mae = mean_absolute_error(y_teste, y_pred)
    rmse = mean_squared_error(y_teste, y_pred, squared = False)
    
    print(modelo,'\n')
    print('Dados de teste')
    print('Acuracia            : {:0.2f}%'.format(acuracia))
    print('MAE                 : {:0.2f}'.format(mae))
    print('MSE                 : {:0.2f}'.format(mse))
    print('RMSE                : {:0.2f}\n'.format(rmse))
In [99]:
def treinaRegressao_GridSearchCV(modelo, params_, x_treino, y_treino, x_teste, y_teste,\
                        n_jobs = -1, cv = 5, refit = True, scoring = None, salvar_resultados = False,\
                        retorna_modelo = False):
    grid = GridSearchCV(modelo, params_, n_jobs = n_jobs, cv = cv, refit = refit, scoring = scoring)
    
    grid.fit(x_treino, y_treino, eval_set = [(x_treino, y_treino)], verbose = False)
    pred = grid.predict(x_teste)
    modelo_ = grid.best_estimator_

    print(grid.best_params_)
    
    reportModeloRegressao(modelo_, x_teste, y_teste) 
    
    if salvar_resultados:
        resultados_df = pd.DataFrame(grid.cv_results_)
        
        if retorna_modelo:
            return resultados_df, modelo_
        else:
            resultados_df
        
    if retorna_modelo:
        return modelo_
In [100]:
# Convertendo dtHistory para serial
dtHistoryIndexado = dtHistory.copy()
dtHistoryIndexado.index = dtHistoryIndexado['Data']
dtHistoryIndexado = dtHistoryIndexado['Vendas']
In [101]:
# Separando base de treino e teste
split_date = '2018-06-16'
train = dtHistoryIndexado.loc[dtHistoryIndexado.index <= split_date].copy()
test = dtHistoryIndexado.loc[dtHistoryIndexado.index > split_date].copy()

# Convertendo base de treio e teste para dataframe
train = pd.DataFrame(train)
test = pd.DataFrame(test)
In [102]:
# Separando X e y em treino e teste
X_train, y_train = create_features(train, label = 'Vendas')
X_test, y_test = create_features(test, label = 'Vendas')
In [103]:
%%time

if EXECUTE_GRID_SEARCH:
    '''params  = {
            'max_depth': np.arange(4, 8, step = 1),
            'gamma': np.arange(1.1, step = 0.1),
            'reg_alpha' : np.arange(1.1, step = 0.1),
            'reg_lambda' : np.arange(1.1, step = 0.1),
            'colsample_bytree' : np.arange(1.1, step = 0.1),
            'min_child_weight' : np.arange(11, step = 1),
            'eta': np.arange(1.1, step = 0.1),
            'n_estimators': [1000],
            'seed': [SEED],
            'early_stopping_rounds': [10],
            'objective': ['reg:squarederror']
        }'''

    params  = {
            'max_depth': np.arange(4, 8, step = 1),
            'gamma': np.arange(1.1, step = 0.1),
            'reg_alpha' : np.arange(1.1, step = 0.1),
            'reg_lambda' : np.arange(1.1, step = 0.1),
            'eta': np.arange(1.1, step = 0.1),
            'n_estimators': [1000],
            'seed': [SEED],
            'early_stopping_rounds': [10],
            'objective': ['reg:squarederror']
        }

    model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', seed = SEED)

    treinaRegressao_GridSearchCV(model, params, X_train, y_train, X_test, y_test, scoring = 'neg_root_mean_squared_error')

{'early_stopping_rounds': 10, 'eta': 0.1, 'gamma': 0.7000000000000001, 'max_depth': 4, 'n_estimators': 1000, 'objective': 'reg:squarederror', 'reg_alpha': 0.7000000000000001, 'reg_lambda': 0.9, 'seed': 84796315}
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
             colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
             early_stopping_rounds=10, enable_categorical=False, eta=0.1,
             eval_metric=None, gamma=0.7000000000000001, gpu_id=-1,
             grow_policy='depthwise', importance_type=None,
             interaction_constraints='', learning_rate=0.100000001, max_bin=256,
             max_cat_to_onehot=4, max_delta_step=0, max_depth=4, max_leaves=0,
             min_child_weight=1, missing=nan, monotone_constraints='()',
             n_estimators=1000, n_jobs=0, num_parallel_tree=1, predictor='auto',
             random_state=84796315, reg_alpha=0.7000000000000001, ...) 

Dados de teste
Acuracia            : 71.67%
MAE                 : 59.34
MSE                 : 4667.55
RMSE                : 68.32

CPU times: total: 1h 45min 56s
Wall time: 1h 45min 37s
In [104]:
# Treinando modelo
model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', max_depth = 4, gamma = 0.7, reg_alpha = 0.7, reg_lambda = 0.9, eta = 0.1, seed = SEED)
model.fit(X_train, y_train,
        eval_set = [(X_train, y_train)],
        verbose = False)
Out[104]:
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
             colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
             early_stopping_rounds=10, enable_categorical=False, eta=0.1,
             eval_metric=None, gamma=0.7, gpu_id=-1, grow_policy='depthwise',
             importance_type=None, interaction_constraints='',
             learning_rate=0.100000001, max_bin=256, max_cat_to_onehot=4,
             max_delta_step=0, max_depth=4, max_leaves=0, min_child_weight=1,
             missing=nan, monotone_constraints='()', n_estimators=1000,
             n_jobs=0, num_parallel_tree=1, predictor='auto',
             random_state=84796315, reg_alpha=0.7, ...)
In [105]:
plot_importance(model, height = 1)
Out[105]:
<AxesSubplot:title={'center':'Feature importance'}, xlabel='F score', ylabel='Features'>
In [106]:
test['Previsões'] = model.predict(X_test)
XGBresults = pd.concat([test, train], sort = False)
In [107]:
pred = test['Previsões']
test_ = test['Vendas']
In [108]:
mse, rmse, mae, mape, mtv, mpv = calculate_metrics_off(test_, pred)
In [109]:
result = record(result, 'Original - XGBoost', mse, rmse, mae, mape, mtv, mpv, False)
In [110]:
# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
                        y = seriesOriginal,
                        name = 'Real'),
             go.Scatter(x = test.index,
                        y = test['Previsões'],
                        name = 'Previsto')]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'Original - XGBoost')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'xgboost_original', auto_open = False)
In [111]:
result
Out[111]:
Algorithm MSE RMSE MAE MAPE Mean_True_Value Mean_Predict_Value
0 Original - Time Series Regression 5296.329792 60.784585 42.960631 25.553663 179.560484 168.008596
1 Deseasonal - Time Series Regression 5081.77776 58.208746 40.122247 23.488563 179.560484 164.670085
2 BoxCox - Time Series Regression 5328.602465 61.060147 43.254815 25.416806 179.560484 165.399796
3 Stationary - Time Series Regression 6057.828087 64.323443 46.538896 27.899678 179.560484 184.267746
4 Original - Exponential Smoothing 4749.165741 55.864659 38.516078 22.095183 179.560484 165.919134
5 Deseasonal - Exponential Smoothing 4774.382729 55.5275 37.488666 21.8423 179.560484 166.749778
6 BoxCox - Exponential Smoothing 4653.246555 54.946508 37.555881 21.63599 179.560484 166.093056
7 Stationary - Exponential Smoothing 5658.627363 63.86568 43.658278 26.080376 179.560484 179.444888
8 Original - Auto Arima 5707.329361 63.881702 45.981605 27.441516 179.560484 169.582966
9 Deseasonal - Auto Arima 9591.527185 77.052424 59.641758 35.25454 179.560484 191.567949
10 BoxCox - Auto Arima 11358.34545 84.301119 65.970324 39.351481 179.560484 206.5619
11 Stationary - Auto Arima 9014.243165 73.755288 56.443311 33.402952 179.560484 192.709473
12 Original - Long Short Term Memory (LSTM) 48171.637513 139.146814 101.418725 70.382288 179.560484 137.450098
13 Deseasonal - Long Short Term Memory (LSTM) 4989.531811 56.707357 38.964999 22.768957 179.560484 165.563544
14 BoxCox - Long Short Term Memory (LSTM) 5341.400913 62.574598 44.292526 26.569442 179.560484 171.38447
15 Stationary - Long Short Term Memory (LSTM) 16091.513901 94.546789 74.40401 43.338943 179.560484 220.839974
16 Original - XGBoost 4667.546621 68.319445 59.339327 28.327163 221.854839 190.8423

4. Comparação¶

In [112]:
result.to_csv('export\\resultados_time_series.csv')
In [113]:
# Tratando nomes e criando colunas de controle para os resultados gerados
topResult = (
    result 
    .assign(Full_Name = lambda x: x.Algorithm.apply(lambda a: a.split('(')[0]
                                                   .replace('ARIMA', 'Auto Arima')
                                                   .replace('Long Short Term Memory', 'LSTM')))
    .assign(Data_Category = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[0]))
    .assign(Algorithm = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[1].split('(')[0]
                                                   .replace('ARIMA', 'Auto Arima')
                                                   .replace('Long Short Term Memory', 'LSTM')))
    .sort_values('MSE')
)
In [114]:
topResult.head()
Out[114]:
Algorithm MSE RMSE MAE MAPE Mean_True_Value Mean_Predict_Value Full_Name Data_Category
6 Exponential Smoothing 4653.246555 54.946508 37.555881 21.63599 179.560484 166.093056 BoxCox - Exponential Smoothing BoxCox
16 XGBoost 4667.546621 68.319445 59.339327 28.327163 221.854839 190.8423 Original - XGBoost Original
4 Exponential Smoothing 4749.165741 55.864659 38.516078 22.095183 179.560484 165.919134 Original - Exponential Smoothing Original
5 Exponential Smoothing 4774.382729 55.5275 37.488666 21.8423 179.560484 166.749778 Deseasonal - Exponential Smoothing Deseasonal
13 LSTM 4989.531811 56.707357 38.964999 22.768957 179.560484 165.563544 Deseasonal - LSTM Deseasonal
In [115]:
# Plot dos resultados obtidos por ordem ascendente do MSE

colors = {'Time Series Regression':'#1688F2',
          'Exponential Smoothing':'#17FC2E',
          'Auto Arima': '#E6C220',
          'LSTM ': '#FC4417',
           'XGBoost': '#A13BF5'}

# plotly figure
fig = go.Figure(layout = go.Layout(yaxis = {'type': 'category', 'title': 'Algoritmo e Categoria'},
                        xaxis = {'title': 'MSE'},
                        title = 'MSE por Algoritmo e Tipo de Dado'))

for t in topResult['Algorithm'].unique():
    topResultFiltered = topResult[topResult['Algorithm']== t].copy()
    fig.add_traces(go.Bar(x = topResultFiltered['MSE'], y = topResultFiltered['Full_Name'], name = str(t),\
                          marker_color = str(colors[t]), orientation = 'h',
                          text = round(topResultFiltered['MSE'].astype(np.double)), textposition = 'outside'))
    
    
fig.update_layout(yaxis = {'categoryorder':'total descending'}, autosize = False,
                  width = 1450,
                  height = 800)
    
fig.show()

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'mse_algoritmo_tipo_dado', auto_open = False)

4. Previsões Futuras¶

In [116]:
pred_range = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 10, 31))
split_range = TimeSeriesSplit(n_splits = 2, max_train_size = pred_range.shape[0], test_size = pred_range.shape[0])

4.1 BoxCox - Exponential Smoothing¶

In [117]:
BATCH_SIZE = 1000
In [118]:
# Alocando melhor modelo a memória e separando base de treino
data = seriesBoxCox.copy()
train = data[datetime.date(2017, 1, 1): datetime.date(2018, 6, 16)]
In [119]:
# Treinando modelo baseado dos parâmetros descobertos na fase de modelagem
ES = (
    ExponentialSmoothing(train, trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
    .fit(smoothing_level = 0.1, smoothing_trend = 0.7, smoothing_seasonal = 0.0, method = 'ls', damping_trend = 0.8)
)
In [120]:
# Calculando a previsão até o final do ano de 2018
pred = ES.predict(str(data.index[0]), '2018-12-31')
pred, data, _ = calcPredTestBoxCox(train, pred, data)
In [121]:
# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
                        y = seriesOriginal,
                        name = 'Real'),
             go.Scatter(x = pred.index,
                        y = pred,
                        name = 'Previsto')]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'BoxCox - Exponential Smoothing')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31', 
              annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'green', opacity = 0.2, line_width = 0)

fig.update_layout(yaxis_range = [0 , 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'exponential_smoothing_boxcox', auto_open = False)
In [122]:
ES.summary()
Out[122]:
ExponentialSmoothing Model Results
Dep. Variable: None No. Observations: 532
Model: ExponentialSmoothing SSE 7800.497
Optimized: True AIC 1452.579
Trend: Additive BIC 1503.899
Seasonal: Additive AICC 1453.392
Seasonal Periods: 7 Date: Tue, 22 Nov 2022
Box-Cox: False Time: 11:40:35
Box-Cox Coeff.: None
coeff code optimized
smoothing_level 0.1000000 alpha False
smoothing_trend 0.7000000 beta False
smoothing_seasonal 0.000000 gamma False
initial_level 3.5625210 l.0 True
initial_trend 1.3021796 b.0 True
damping_trend 0.8000000 phi False
initial_seasons.0 -8.0164149 s.0 True
initial_seasons.1 -3.7667775 s.1 True
initial_seasons.2 -3.8599627 s.2 True
initial_seasons.3 -4.3632835 s.3 True
initial_seasons.4 -5.0029325 s.4 True
initial_seasons.5 -6.0412289 s.5 True
initial_seasons.6 -9.3733856 s.6 True

4.2 Deseasonal - LSTM¶

In [123]:
BATCH_SIZE = 30
In [124]:
# Definindo dataset
dtHistoryLSTM = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date

# Separando periodo minimo e maximo
startLSTM = dtHistoryLSTM.min()
endLSTM = dtHistoryLSTM.max()

# Separando IDs das lacunas
idxLSTM = pd.date_range(startLSTM, endLSTM, normalize = True)

# Transformando dt para série e realizando contagem de valores diários
seriesOriginalLSTM = dtHistoryLSTM.value_counts(sort = False).sort_index().reindex(idxLSTM, fill_value = 0)

# Removendo outliers
seriesOriginalLSTM = seriesOriginalLSTM[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]

# Adicionando  predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
seriesOriginalTSTRPredictions = seriesOriginalLSTM.append(newPredictions)

# Gerando STL para separar o Deseasonal
stl = STL(seriesOriginalTSTRPredictions)
res = stl.fit()

# Separando Deseasonal
seriesOriginalTSTRPredictionsDeseasonal = res.observed - res.seasonal

# Separando index da série
dataTime = seriesOriginalTSTRPredictionsDeseasonal.copy().index

# Separando base de treino e teste
train_size = 532
test_size = len(seriesOriginalTSTRPredictionsDeseasonal) - train_size
train = seriesOriginalTSTRPredictionsDeseasonal[0:train_size]
test = seriesOriginalTSTRPredictionsDeseasonal[train_size:len(seriesOriginalTSTRPredictionsDeseasonal)]
In [125]:
# Normalização e reshape do conjunto de treino
scaler = MinMaxScaler(feature_range=(0, 1))
scaler.fit(train.values.reshape(-1, 1))
train_input = scaler.transform(train.values.reshape(-1, 1))

# Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
test_input = train_input[-(FEATURES + 1):]
test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)

# Reset da sessão
clear_session()

# Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
initializer = GlorotUniform(seed = SEED)
model = Sequential()

# 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
model.add(Dense(1, kernel_initializer = initializer))
model.add(BatchNormalization())
model.add(Dense(1, kernel_initializer = initializer))

# Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')

# Treinando o modelo com otimizador Adam
model.compile(loss = 'mse', optimizer = Adam(learning_rate = 0.01), metrics = ['mae'])
model.fit(train_gen, epochs = 100, callbacks = [early_stopping], verbose = True)
pred = np.empty(test.shape[0])

# Realizando predições no conjunto de teste
for i in range(len(pred)):
    prediction = model.predict(test_gen, verbose = False)
    pred[i] = prediction
    #print(test_gen[0][0], ' = ', prediction)
    test_input = np.delete(test_input, 0, 0)
    test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
    test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)

# Reorganizando o shape e chamando a função de calculo
pred = scaler.inverse_transform(pred.reshape(-1,1)).reshape(-1)
pred, test, _ = calcPredTestseriesDeseasonal(train, pred, test)

Epoch 1/100
18/18 [==============================] - 2s 4ms/step - loss: 0.0231 - mae: 0.0996
Epoch 2/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0083 - mae: 0.0699
Epoch 3/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0067 - mae: 0.0649
Epoch 4/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0071 - mae: 0.0647
Epoch 5/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0069 - mae: 0.0666
Epoch 6/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0066 - mae: 0.0637
Epoch 7/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0071 - mae: 0.0652
Epoch 8/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0068 - mae: 0.0670
Epoch 9/100
18/18 [==============================] - 0s 4ms/step - loss: 0.0067 - mae: 0.0660
In [126]:
# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = dataTime,
                        y = seriesOriginalTSTRPredictionsDeseasonal[:594],
                        name = 'Real'),
             go.Scatter(x = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 12, 31)),
                        y = pred,
                        name = 'Previsto')]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'Deseasonal - LSTM')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31', 
              annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'green', opacity = 0.2, line_width = 0)

fig.update_layout(yaxis_range = [0 , 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'lstm_deseasonal', auto_open = False)

4.3 Deseasonal - Time Series Regressor¶

In [127]:
BATCH_SIZE = 1000
In [128]:
# Definindo dataset
dtHistoryTSR = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date

# Separando periodo minimo e maximo
startTSR = dtHistoryTSR.min()
endTSR = dtHistoryTSR.max()

# Separando IDs das lacunas
idxTSR = pd.date_range(startTSR, endTSR, normalize = True)

# Transformando dt para série e realizando contagem de valores diários
seriesOriginalTSR = dtHistoryTSR.value_counts(sort = False).sort_index().reindex(idxTSR, fill_value = 0)

# Removendo outliers
seriesOriginalTSR = seriesOriginalTSR[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]

# Adicionando  predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
seriesOriginalTSRRPredictions = seriesOriginalTSR.append(newPredictions)

# Gerando STL para separar o Deseasonal
stl = STL(seriesOriginalTSRRPredictions)
res = stl.fit()

# Separando Deseasonal
seriesOriginalTSRRPredictionsDeseasonal = res.observed - res.seasonal

# Separando base de treino e teste
train_size = 532
test_size = len(seriesOriginalTSRRPredictionsDeseasonal) - train_size
train = seriesOriginalTSRRPredictionsDeseasonal[0:train_size]
test = seriesOriginalTSRRPredictionsDeseasonal[train_size:len(seriesOriginalTSRRPredictionsDeseasonal)]
In [129]:
# Gerando generator de treino
gen = TimeseriesGenerator(train, train, FEATURES, batch_size = BATCH_SIZE)

X_train = gen[0][0]
y_train = gen[0][1]

# Treinando modelo
lr = LinearRegression()
lr.fit(X_train, y_train)

# Separando primeira leva de treino
X_pred = y_train[-FEATURES:].reshape(1,-1)
pred = np.empty(test.shape[0])

# Realizando predições e realocando vetor de entrada
for i in range(len(pred)):
    forecast = lr.predict(X_pred)
    X_pred = np.delete(X_pred, 0, 1)
    X_pred = np.concatenate((X_pred, forecast.reshape(-1, 1)), 1)    
    pred[i] = forecast

# Convertendo deseasonal de volta para os valores originais
pred, test, _ = calcPredTestseriesDeseasonal(train, pred, test)
In [130]:
# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = dataTime,
                        y = seriesOriginalTSTRPredictionsDeseasonal[:594],
                        name = 'Real'),
             go.Scatter(x = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 12, 31)),
                        y = pred,
                        name = 'Previsto')]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'Deseasonal - Time Series Regressor')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31', 
              annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'green', opacity = 0.2, line_width = 0)

fig.update_layout(yaxis_range = [0 , 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'time_series_regressor_deseasonal', auto_open = False)

4.4 BoxCox - Arima¶

In [131]:
# Definindo dataset
dtHistoryArima = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date

# Separando periodo minimo e maximo
startArima = dtHistoryArima.min()
endArima = dtHistoryArima.max()

# Separando IDs das lacunas
idxArima = pd.date_range(startArima, endArima, normalize = True)

# Transformando dt para série e realizando contagem de valores diários
seriesOriginalArima = dtHistoryArima.value_counts(sort = False).sort_index().reindex(idxArima, fill_value = 0)

# Removendo outliers
seriesOriginalArima = seriesOriginalArima[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]

# Adicionando  predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
seriesOriginalArimaPredictions = seriesOriginalArima.append(newPredictions)

# Gerando STL para separar cada um dos componentes
stl = STL(seriesOriginalArimaPredictions)
res = stl.fit()

# Separando boxcox
seriesBoxCoxArima, lmbda = boxcox(seriesOriginalArimaPredictions + 1)
seriesBoxCoxArima = pd.Series(seriesBoxCoxArima, index = seriesOriginalArimaPredictions.index)

# Separando base de treino e teste
train_size = 532
test_size = len(seriesBoxCoxArima) - train_size
train = seriesBoxCoxArima[0:train_size]
test = seriesBoxCoxArima[train_size:len(seriesBoxCoxArima)]
In [132]:
arimaModel = auto_arima(seriesBoxCoxArima, m = FEATURES, seasonal = True)
arimaModel
Out[132]:
ARIMA(order=(1, 1, 1), scoring_args={}, seasonal_order=(1, 0, 2, 7),
      suppress_warnings=True, with_intercept=False)
In [133]:
arimaModel.fit(train)
pred = arimaModel.predict(test.shape[0])

pred, test, _ = calcPredTestBoxCox(train, pred, test)
In [134]:
# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = pd.date_range(datetime.date(2017, 1, 1), datetime.date(2018, 8, 17)),
                        y = seriesOriginalArima[:594],
                        name = 'Real'),
             go.Scatter(x = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 12, 31)),
                        y = pred,
                        name = 'Previsto')]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'BoxCox - Arima')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31', 
              annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'green', opacity = 0.2, line_width = 0)

fig.update_layout(yaxis_range = [0 , 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'arima_boxcox', auto_open = False)

4.5 Original - XGBoost¶

In [135]:
# Convertendo dtHistory para serial
dtHistoryIndexado = dtHistory.copy()
dtHistoryIndexado.index = dtHistoryIndexado['Data']
dtHistoryIndexado = dtHistoryIndexado['Vendas']
In [136]:
# Separando base de treino e teste
split_date = '2018-06-16'
train = dtHistoryIndexado.loc[dtHistoryIndexado.index <= split_date].copy()
test = dtHistoryIndexado.loc[dtHistoryIndexado.index > split_date].copy()

# Convertendo base de treio e teste para dataframe
train = pd.DataFrame(train)
test = pd.DataFrame(test)
In [137]:
# Separando X e y em treino e teste
X_train, y_train = create_features(train, label = 'Vendas')
X_test, y_test = create_features(test, label = 'Vendas')
In [138]:
# Treinando modelo
model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', max_depth = 4, gamma = 0.7, reg_alpha = 0.7, reg_lambda = 0.9, eta = 0.1, seed = SEED)
model.fit(X_train, y_train,
        eval_set = [(X_train, y_train)],
        verbose = False)
Out[138]:
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
             colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
             early_stopping_rounds=10, enable_categorical=False, eta=0.1,
             eval_metric=None, gamma=0.7, gpu_id=-1, grow_policy='depthwise',
             importance_type=None, interaction_constraints='',
             learning_rate=0.100000001, max_bin=256, max_cat_to_onehot=4,
             max_delta_step=0, max_depth=4, max_leaves=0, min_child_weight=1,
             missing=nan, monotone_constraints='()', n_estimators=1000,
             n_jobs=0, num_parallel_tree=1, predictor='auto',
             random_state=84796315, reg_alpha=0.7, ...)
In [139]:
# Adicionando  predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
newPredictions = pd.DataFrame(newPredictions)
newPredictions = create_features(newPredictions, label = None)
newPredictionsTest = pd.concat([X_test, newPredictions])
test = pd.concat([test, newPredictions])
In [140]:
test['Previsões'] = model.predict(newPredictionsTest)
In [141]:
# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
                        y = seriesOriginal,
                        name = 'Real'),
             go.Scatter(x = test.index,
                        y = test['Previsões'],
                        name = 'Previsto')]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'Original - XGBoost')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31', 
              annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'green', opacity = 0.2, line_width = 0)

fig.update_layout(yaxis_range = [0 , 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'original_xgboost_projecao', auto_open = False)

5. Melhor Modelo¶

In [142]:
# Convertendo dtHistory para serial
dtHistoryIndexado = dtHistory.copy()
dtHistoryIndexado.index = dtHistoryIndexado['Data']
dtHistoryIndexado = dtHistoryIndexado['Vendas']
In [143]:
# Separando base de treino e teste
split_date = '2018-06-16'
train = dtHistoryIndexado.loc[dtHistoryIndexado.index <= split_date].copy()
test = dtHistoryIndexado.loc[dtHistoryIndexado.index > split_date].copy()

# Convertendo base de treio e teste para dataframe
train = pd.DataFrame(train)
test = pd.DataFrame(test)
In [144]:
# Separando X e y em treino e teste
X_train, y_train = create_features(train, label = 'Vendas')
X_test, y_test = create_features(test, label = 'Vendas')
In [145]:
# Treinando modelo
model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', max_depth = 4, gamma = 0.7, reg_alpha = 0.7, reg_lambda = 0.9, eta = 0.1, seed = SEED)
model.fit(X_train, y_train,
        eval_set = [(X_train, y_train)],
        verbose = False)
Out[145]:
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
             colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
             early_stopping_rounds=10, enable_categorical=False, eta=0.1,
             eval_metric=None, gamma=0.7, gpu_id=-1, grow_policy='depthwise',
             importance_type=None, interaction_constraints='',
             learning_rate=0.100000001, max_bin=256, max_cat_to_onehot=4,
             max_delta_step=0, max_depth=4, max_leaves=0, min_child_weight=1,
             missing=nan, monotone_constraints='()', n_estimators=1000,
             n_jobs=0, num_parallel_tree=1, predictor='auto',
             random_state=84796315, reg_alpha=0.7, ...)
In [146]:
newPredictionsTest = pd.concat([X_train, X_test])
test = pd.concat([train, test])
In [147]:
test['Previsões'] = model.predict(newPredictionsTest)
In [148]:
test['ci'] = 0.04 * test['Vendas']
test['ci_lower'] = test['Vendas'] - test['ci']
test['ci_upper'] = test['Vendas'] + test['ci']
In [149]:
test.tail()
Out[149]:
Vendas date hour dayofweek quarter month year dayofyear dayofmonth weekofyear Previsões ci ci_lower ci_upper
Data
2018-08-13 284 2018-08-13 0 0 3 8 2018 225 13 33 239.025238 11.36 272.64 295.36
2018-08-14 311 2018-08-14 0 1 3 8 2018 226 14 33 231.174133 12.44 298.56 323.44
2018-08-15 282 2018-08-15 0 2 3 8 2018 227 15 33 255.116379 11.28 270.72 293.28
2018-08-16 316 2018-08-16 0 3 3 8 2018 228 16 33 237.553925 12.64 303.36 328.64
2018-08-17 249 2018-08-17 0 4 3 8 2018 229 17 33 203.458313 9.96 239.04 258.96
In [150]:
# Plot

# Definição dos dados no plot
plot_data = [go.Scatter(x = test.index,
             y = test['Previsões'],
             name = 'Previsto'),
            go.Scatter(
                name = '95% CI Upper',
                x = test.index,
                y = test['ci_upper'],
                mode = 'lines',
                marker = dict(color = '#444'),
                line = dict(width = 0),
                showlegend = False
            ),
            go.Scatter(
                name = '95% CI Lower',
                x = test.index,
                y = test['ci_upper'],
                marker = dict(color = '#444'),
                line = dict(width = 0),
                mode = 'lines',
                fillcolor = 'rgba(68, 68, 68, 0.3)',
                fill = 'tonexty',
                showlegend = False
            )
        ]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
                        yaxis = {'title': 'Vendas'}, 
                        title = 'Original - XGBoost')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

fig.update_layout(yaxis_range = [0 , 800])

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'xgboost_original_base_teste', auto_open = False)
In [151]:
fig = go.Figure([
    go.Scatter(
        name = 'Previsto',
        x = test.index,
        y = test['Previsões'],
    ),
    go.Scatter(
        name = '95% CI Upper',
        x = test.index,
        y = test['ci_upper'],
        mode = 'lines',
        marker = dict(color = '#444'),
        line = dict(width = 0),
        showlegend = True
    ),
    go.Scatter(
        name = '95% CI Lower',
        x = test.index,
        y = test['ci_lower'],
        marker = dict(color = '#444'),
        line = dict(width = 0),
        mode = 'lines',
        fillcolor = 'rgba(255, 0, 0, 0.80)',
        fill = 'tonexty',
        showlegend = True
    )
])

fig.update_layout(
    xaxis_title = 'Período',
    yaxis_title = 'Vendas',
    title = 'Original - XGBoost',
    hovermode = 'x',
    yaxis_range = [0 , 800]
)

fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17', 
              annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
              annotation = dict(font_size = 23, font_family = 'Times New Roman'),
              fillcolor = 'orange', opacity = 0.2, line_width = 0)

fig.update_layout(yaxis_range = [0 , 800])

fig.update_yaxes(rangemode='tozero')

pyoff.iplot(fig)

# Upando no dash
if upload_to_dash:
    py.plot(fig, filename = 'xgboost_original_limite_confianca', auto_open = False)